C  MODPATH release: Version 4.00 (V4, Release 1, 2-2000)
C    Changes to work with MODFLOW-2000
C
C MODPATH Version 3.00 (V3, Release 2, 5-99)
C Changes:
C   No change from previous release: (V3, Release 1, 9-94)
C***** SUBROUTINES *****
C     DATIN
C***********************
 
C***** SUBROUTINE *****
      SUBROUTINE DATIN (LAYCON,NCON,HEAD,DELT,XMAX,YMAX,DELR,DELC,
     1ZTOP,ZBOT,IBOUND,POR,QX,QY,QZ,IUNIT,IRP,NPRT,FRAC,BUFF,
     2IREV,IBUFF,JLC,ILC,KLC,XLC,YLC,ZLLC,ITMS,QSS,QSTO,ISNK,ITREAD,
     3IZSTOP,IPEP,MAXSTP,
     4INILOC,TSTOP,TIMREL,KKPER,KKSTP,PERLEN,NUMTS,
     5TIMX,IBSTRT,ZLC,HDRY,HNOFLO,TOT,ICMPCT,TBGABS,LAYCBD)
C
      INCLUDE 'idat1.inc'
      DIMENSION LAYCBD(200)
      CHARACTER*256 FLPART,FIL,OUT1,OUT3,FNAME,FNAMS  ! emrl 80 to 256
      CHARACTER*78 MES,STRING,CFACES
      CHARACTER*24 ANAME
      CHARACTER*132 LINE
      CHARACTER*81 LINE2
      COMMON /IDAT2/ IRCHTP,IEVTTP,NRCHOP,NEVTOP
      COMMON /DFNAMS/ FNAMS(6)
C
      DIMENSION RARRAY(10),IARRAY(10)
C
      DIMENSION LAYCON(NLAY),NCON(NLAY),DELR(NCOL),DELC(NROW),
     1XMAX(NCOL),YMAX(NROW),ZTOP(NZDIM),ZBOT(NZDIM),
     2HEAD(NCOL,NROW,NLAY),IBOUND(NCOL,NROW,NLAY),POR(NCOL,NROW,NLPOR),
     3QX(NCP1,NROW,NLAY),QY(NCOL,NRP1,NLAY),QZ(NCOL,NROW,NLP1),
     4IUNIT(NUNIT),BUFF(NCOL,NROW,NLAY),IBUFF(NCOL,NROW,NLAY),
     5JLC(NPART),ILC(NPART),KLC(NPART),XLC(NPART),YLC(NPART),
     6ZLLC(NPART),QSS(NCOL,NROW,NLAY),
     8QSTO(NCOL,NROW,NLAY),INILOC(NPART),
     9PERLEN(NPER),NUMTS(NPER),TIMX(NPER),
     1IBSTRT(NCOL,NROW,NLAY),TOT(NPART),ZLC(NPART)
C
C
      ITREAD=0
C
C  Read DIS file
      CALL DIS2RP(IUNIT(1),I7,DELR,DELC,ZTOP,ZBOT,LAYCBD,NCON,
     1          PERLEN,NUMTS,TIMX,ISS,NCOL,NROW,NLAY,NPER)
      XMAX(1)=DELR(1)
      IF(NCOL.GT.1) THEN
         DO 40 J=2,NCOL
         XMAX(J)=XMAX(J-1)+DELR(J)
40       CONTINUE
      END IF
C
      YMAX(NROW)=DELC(NROW)
      IF(NROW.GT.1) THEN
         DO 50 I=2,NROW
         II=NROW+1-I
         YMAX(II)=YMAX(II+1)+DELC(II)
50       CONTINUE
      END IF
C
      IU=I1
      NPRT=NPART
C
C  WRITE NCOL, NROW, NLAY, NCBL, AND IGRID TO SUMMARY.PTH
C
      WRITE(I7,5000)
5000  FORMAT('  ')
      WRITE(I7,5010) NCOL,NROW,NLAY,NCBL
5010  FORMAT(I5,' COLUMNS',I5,' ROWS',I5,' LAYERS',I5,' CONFINING LAYERS
     1')
      WRITE(I7,5020) IGRID
5020  FORMAT(' IGRID (GRID TYPE CODE) IS',I2)
C
C  STEADY OR TRANSIENT?
C
C... IF STEADY STATE THEN...
      IF(ISS.EQ.1) THEN
      NPER=1
      INHQ=1
      KKPER=1
      KKSTP=1
      TBGABS=0.0
      TBEGIN=0.0
      MES= 'DO YOU WANT TO STOP COMPUTING PATHS AFTER A SPECIFIED LENGTH
     1 OF TIME ?'
      CALL YESNO(MES,IANS,'2.1.41')
      TSTOP= 1.0E+30
      IF(IANS.EQ.1) THEN
      CALL MPRMPT
     1('ENTER: MAXIMUM TRACKING TIME & TIME UNITS CONVERSION FACTOR')
      CALL GETR(2,RARRAY,0,0.,0.,'2.1.30')
      TSTOP= RARRAY(1)*RARRAY(2)
      END IF
 
c... else if transient then...
      ELSE
      IF(NPER.EQ.0) THEN
      WRITE(*,*)
     1'NUMBER OF STRESS PERIODS (NPER) SET = 0 IN MAIN DATA FILE. STOP.'
		call stopfile  ! emrl
      STOP
      END IF
      CALL MPRMPT('DEFINE A REFERENCE TIME FOR RELEASING PARTICLES ...')
      CALL MPRMPT('  SELECT AN OPTION:')
      CALL MPRMPT
     1('     1 = SPECIFY BY ENTERING A STRESS PERIOD AND TIME STEP')
      CALL MPRMPT
     1('     2 = SPECIFY BY ENTERING A VALUE OF SIMULATION TIME')
      CALL GETI(1,IANS,3,1,2,'2.1.2A')
C
      IF(IANS.EQ.2) THEN
      KKPER=0
      KKSTP=0
      CALL MPRMPT
     1('  ENTER: REFERENCE TIME  &  TIME UNITS CONVERSION FACTOR')
      CALL GETR(2,RARRAY,0,0.,0.,'2.1.31')
      TBGABS= RARRAY(1)*RARRAY(2)
      ELSE
      CALL MPRMPT
     1('  ENTER: STRESS PERIOD & TIME STEP ')
      CALL GETI(2,IARRAY,0,0,0,'2.1.2')
      KKPER=IARRAY(1)
      KKSTP=IARRAY(2)
      CALL MPRMPT
     1('  ENTER: RELATIVE TIME WITHIN TIME STEP')
      CALL MPRMPT
     1('         (VALUE FROM 0 TO 1)')
      CALL GETR(1,TIMREL,3,0.,1.,'2.1.32')
      END IF
C
      MES='STOP COMPUTING PATHS AT A SPECIFIED VALUE OF TRACKING TIME ?'
      CALL YESNO(MES,IANS,'2.1.42')
      TSTOP= 1.0E+30
      IF(IANS.EQ.1) THEN
      CALL MPRMPT
     1('ENTER: MAXIMUM TRACKING TIME & TIME UNITS CONVERSION FACTOR')
      CALL GETR(2,RARRAY,0,0.,0.,'2.1.33')
      TSTOP= RARRAY(1)*RARRAY(2)
      END IF
      CALL MPRMPT
     1('SPECIFY AN OPTION FOR READING HEAD AND FLOW RATE DATA:')
      CALL MPRMPT
     1('  1 = READ STANDARD MODFLOW UNFORMATTED FILES & GENERATE A')
      CALL MPRMPT
     1('      COMPOSITE BUDGET FILE')
      CALL MPRMPT
     1('  2 = READ FROM AN EXISTING COMPOSITE BUDGET FILE')
      CALL GETI(1,IANS,3,1,2,'2.1.3')
      IF(IANS .EQ. 1) THEN
        INHQ=3
      ELSE IF(IANS .EQ. 2) THEN
        INHQ=2
      END IF
      IF(FNAMS(1).EQ.' ') THEN
      CALL MPRMPT('ENTER A NAME FOR THE COMPOSITE BUDGET FILE (CBF):')
      CALL GETSTR(FNAME,0,'2.1.36')
      ELSE
        FNAME=FNAMS(1)
      END IF
 
      END IF
C
C  ENTER MODE DATA
C
      CALL MPRMPT('SELECT THE OUTPUT MODE:')
      CALL MPRMPT('    1 = ENDPOINTS')
      CALL MPRMPT('    2 = PATHLINE')
      CALL MPRMPT('    3 = TIME SERIES')
      CALL GETI(1,MODE,3,1,3,'2.1.6')
      MODE=MODE-1
C
C  ENTER TIME STEP INFORMATION FOR MODES 1 AND 2, AND DUMMY VALUES
C  FOR MODE 0
C
      IF (MODE.EQ.0) THEN
      ITMS=0
      END IF
      IF(MODE.EQ.1.OR.MODE.EQ.2) THEN
 
      IF (MODE.EQ.1) THEN
      MES= 'DO YOU WANT TO COMPUTE LOCATIONS AT SPECIFIC POINTS IN TIME
     1?'
      CALL YESNO (MES,ITMS,'2.1.43')
      ELSE IF (MODE.EQ.2) THEN
      ITMS=1
      END IF
 
      IF(ITMS.EQ.1) THEN
      CALL MPRMPT('HOW SHOULD POINTS IN TIME BE SPECIFIED ?')
      CALL MPRMPT('    1 = WITH A CONSTANT TIME INTERVAL')
      CALL MPRMPT('    2 = VALUES OF TIME POINTS ARE READ FROM A FILE')
      CALL GETI(1,ITREAD,3,1,2,'2.1.7')
      ITREAD=ITREAD-1
 
      IF (ITREAD.EQ.1) THEN
        IF(FNAMS(3).EQ.' ') THEN
          CALL MPRMPT('ENTER THE NAME OF THE DATA FILE:')
          CALL GETSTR(FIL,-1,'2.1.38')
          CALL OPNFIL (I8,FIL,1,I7,IBATCH,1)
          CALL PUTSTR(FIL)
        ELSE
          CALL OPNFIL (I8,FNAMS(3),1,I7,IBATCH,1)
        END IF
      ELSE
      CALL MPRMPT
     1('ENTER: TIME INTERVAL & TIME UNITS CONVERSION FACTOR')
      CALL GETR(2,RARRAY,0,0.,0.,'2.1.34')
      DELT=RARRAY(1)*RARRAY(2)
      CALL MPRMPT('ENTER THE MAXIMUM NUMBER OF TIME POINTS ALLOWED')
      CALL GETI(1,MAXSTP,1,0,0,'2.1.8')
      END IF
      END IF
 
      END IF
      CALL MPRMPT('HOW ARE STARTING LOCATIONS TO BE ENTERED?')
      CALL MPRMPT('    1 = FROM AN EXISTING DATA FILE')
      CALL MPRMPT
     1('    2 = ARRAYS OF PARTICLES WILL BE GENERATED INTERNALLY')
      CALL GETI(1,IRP,3,1,2,'2.1.9')
C
C  OPEN FILE FOR PARTICLE STARTING LOCATIONS
C
      IF(IRP.EQ.1) THEN
      IF(FNAMS(2).EQ.' ') THEN
        CALL MPRMPT
     1   ('ENTER NAME OF DATA FILE CONTAINING STARTING LOCATIONS:')
        CALL GETSTR(FLPART,-1,'2.1.39')
        CALL OPNFIL (I6,FLPART,1,I7,IBATCH,1)
        CALL PUTSTR(FLPART)
      ELSE
        CALL OPNFIL (I6,FNAMS(2),1,I7,IBATCH,1)
      END IF
 
      ELSE IF (IRP.EQ.2) THEN
      MES= 'DO YOU WANT TO STORE INTERNALLY-GENERATED STARTING LOCATIONS
     1 ON DISK ?'
      CALL YESNO(MES,ILSTOR,'2.1.44')
      IF(ILSTOR.EQ.1) THEN
      CALL MPRMPT('ENTER A FILE NAME:')
      CALL GETSTR(FLPART,-1,'2.1.40')
      CALL OPNFIL (I6,FLPART,2,I7,IBATCH,3)
      CALL PUTSTR(FLPART)
      END IF
      END IF
 
C
C  CHOOSE BETWEEN FORWARD AND REVERSE TRACKING
C
      CALL MPRMPT('IN WHICH DIRECTION SHOULD PARTICLES BE TRACKED?')
      CALL MPRMPT('    1 = FORWARD IN THE DIRECTION OF FLOW')
      CALL MPRMPT('    2 = BACKWARDS TOWARD RECHARGE LOCATIONS')
      CALL GETI(1,IREV,3,1,2,'2.1.10')
      IREV=IREV-1
C
C  SELECT DISCHARGE CRITERION
C
      CALL MPRMPT
     1('HOW SHOULD PARTICLES BE TREATED WHEN THEY ENTER CELLS WITH INTER
     2NAL SINKS ?')
      CALL MPRMPT
     1('    1 = PASS THROUGH WEAK SINK CELLS')
      CALL MPRMPT
     1('    2 = STOP AT WEAK SINK CELLS')
      CALL MPRMPT
     1('    3 = STOP AT WEAK SINK CELLS THAT EXCEED A SPECIFIED STRENGTH
     2')
      CALL GETI(1,ISNK,3,1,3,'2.1.11')
      ISNK=ISNK-1
 
      FRAC=0.99
      IF(ISNK.EQ.2) THEN
      CALL MPRMPT('ENTER A NUMBER BETWEEN 0 AND 1:')
      CALL MPRMPT('    (0.0 => NONE OF THE INFLOW TO THE CELL IS DISCHAR
     1GED TO INTERNAL SINKS)')
      CALL MPRMPT('    (1.0 => ALL INFLOW TO THE CELL IS DISCHARGED TO I
     1NTERNAL SINKS)')
      CALL GETR(1,FRAC,3,0.,1.,'2.1.35')
      IF(FRAC.GE.0.99) FRAC=0.99
      END IF
C
C  SPECIFY AUTOMATIC DISCHARGE ZONE
C
      IZSTOP=0
      IPEP=0
      MES= 'DO YOU WANT TO STOP PARTICLES WHENEVER THEY ENTER ONE SPECIF
     1IC ZONE ?'
      CALL YESNO (MES,IANS,'2.1.45')
      IF(IANS.EQ.1) THEN
      CALL MPRMPT('ENTER THE ZONE NUMBER (MUST BE > 1):')
      CALL GETI(1,IZSTOP,1,2,0,'2.1.12')
 
      IF(MODE.EQ.0) THEN
      CALL MPRMPT('SPECIFY WHICH ENDPOINTS TO RECORD:')
      CALL MPRMPT('   1 = ENDPOINT DATA RECORDED FOR ALL PARTICLES')
      CALL MPRMPT('   2 = ENDPOINT DATA RECORDED ONLY FOR PARTICLES')
      WRITE (STRING,5060) IZSTOP
      CALL MPRMPT(STRING)
5060  FORMAT ('        TERMINATING IN ZONE ',I3)
      CALL GETI(1,IPEP,3,1,2,'2.1.13')
      IPEP=IPEP-1
      ELSE
      IPEP=0
      END IF
 
      END IF
C
      IF (IBATCH.EQ.0)
     1WRITE (*,*) 'READING AND PROCESSING INPUT DATA...'
 
C
C  FIND OUT IF IBOUND ARRAY IS READ IN AS A CROSS SECTION
C
      IXREAD=0
      ICMPCT=0
      READ(IU,'(A)') LINE
      CALL UPCASE(LINE)
      IF(INDEX(LINE,'XSECTION').NE.0) IXREAD=1
      IF(INDEX(LINE,'COMPACT').NE.0) ICMPCT=1
      IF(INDEX(LINE,'BINARY') .NE. 0) ICMPCT=2
      IF(IXREAD.EQ.1 .AND. NROW.GT.1) THEN
        WRITE(I7,*)
     1'<XSECTION> WAS SPECIFIED FOR A MODEL WITH MORE THAN 1 ROW. STOP.'
		call stopfile  ! emrl
        STOP
      END IF
 
C
C  CREATE AND OPEN 'ENDPOINT', 'PATHLINE' AND (OR) 'TIMESERS' FILES
C
 
      IF(ICMPCT.LT.2) THEN
        I3TEMP= I3
        OUT1= 'endpoint'
      ELSE IF(ICMPCT.EQ.2) THEN
        I3TEMP= -I3
        OUT1= 'endpoint.bin'
      END IF
      IF(FNAMS(4).NE.' ') OUT1=FNAMS(4)
      CALL OPNFIL (I3TEMP,OUT1,4,I7,IBATCH,3)
 
      IF(MODE.EQ.1) THEN
        IF(ICMPCT.LT.2) THEN
          I2TEMP= I2
          OUT1= 'pathline'
        ELSE IF(ICMPCT.EQ.2) THEN
          I2TEMP= -I2
          OUT1= 'pathline.bin'
        END IF
        IF(FNAMS(5).NE.' ') OUT1=FNAMS(5)
        CALL OPNFIL (I2TEMP,OUT1,4,I7,IBATCH,3)
 
      ELSE IF(MODE.EQ.2) THEN
        IF(ICMPCT.LT.2) THEN
          I4TEMP= I4
          OUT3= 'timesers'
        ELSE IF(ICMPCT.EQ.2) THEN
          I4TEMP= -I4
          OUT3= 'timesers.bin'
        END IF
        IF(FNAMS(6).NE.' ') OUT3=FNAMS(6)
        CALL OPNFIL (I4TEMP,OUT3,4,I7,IBATCH,3)
      END IF
 
C
C  LAYER TYPE CODES
C
      CALL UREADI(LAYCON,NLAY,IU,I7)
      WRITE(I7,5000)
      WRITE(I7,*) 'LAYCON (LAYER TYPE CODES):'
      WRITE(I7,5080) (LAYCON(N),N=1,NLAY)
5080  FORMAT(25I3)
      DO 5082 K=1,NLAY
      IF(LAYCON(K).NE.0) LAYCON(K)=3
5082  CONTINUE
C
C
C  IBOUND DATA
C
      ANAME= 'IBOUND'
      IF(IXREAD.EQ.0) THEN
        DO 90 K=1,NLAY
        CALL U2DINT(IBOUND(1,1,K),ANAME,NROW,NCOL,K,IU,I7)
90      CONTINUE
      ELSE
        CALL U2DINT(IBOUND,ANAME,NLAY,NCOL,-1,IU,I7)
      END IF
 
C--- USE IBOUND TO SET THE IBSTRT ARRAY.
      DO 91 K=1,NLAY
      DO 91 I=1,NROW
      DO 91 J=1,NCOL
      IF(ABS(IBOUND(J,I,K)).GT.999) THEN
      WRITE(I7,*)
     1'MODPATH requires the absolute value of IBOUND to be < 1000'
      WRITE(I7,*) 'Stop.'
		call stopfile  ! emrl
      STOP
      END IF
      IF(IBOUND(J,I,K).LT.0) THEN
      IBSTRT(J,I,K) = -IBOUND(J,I,K)
      IBOUND(J,I,K) = -1
      ELSE IF(IBOUND(J,I,K).GT.0) THEN
      IBSTRT(J,I,K)= IBOUND(J,I,K)
      IBOUND(J,I,K) = 1
      ELSE
        IF(LAYCON(K).GT.0) THEN
        IBSTRT(J,I,K)=1
        ELSE
        IBSTRT(J,I,K)=0
        END IF
      END IF
91    CONTINUE
 
C
C  POROSITY DATA
C
      DO 100 K=1,NLAY
      ANAME= 'POROSITY'
      CALL U2DREL(POR(1,1,K),ANAME,NROW,NCOL,K,IU,I7)
      LCON=NCON(K)
      IF(LCON.GT.0) THEN
      KCB=NLAY+LCON
      ANAME= 'CONFINING BED POROSITY'
      CALL U2DREL(POR(1,1,KCB),ANAME,NROW,NCOL,K,IU,I7)
      END IF
100    CONTINUE
C
C  STRESS PERIOD DATA FOR TRANSIENT SIMULATIONS
C  Main data file item 10A.
      IF(ISS.EQ.0) THEN
      READ(IU,'(A)') LINE
      ICOL=1
      CALL URWORD(LINE,ICOL,IWSTRT,IWLAST,3,NDUMMY,TBEGIN,I7,IU)
      WRITE(I7,5320) TBEGIN
5320  FORMAT(1X/1X,'TIME AT BEGINNING OF SIMULATION = ',E14.6)
C
      WRITE(I7,*) 'TIME STEP SUMMARY:'
      TSIM=TBEGIN
      DO 111 N=1,NPER
      KKK=NUMTS(N)
      DO 111 KK=1,KKK
      CALL STPSIZ(PERLEN(N),KK,KKK,TIMX(N),STPL)
      TSIM=TSIM+STPL
      WRITE(I7,5345) N,KK,STPL,TSIM
5345  FORMAT(1X,'PERIOD',I4,' STEP',I4,': DELT =',E14.6,
     1' TIME AT END OF STEP =',E14.6)
111   CONTINUE
      WRITE(I7,*)
C
      IF(KKPER.GT.0 .AND. KKSTP.GT.0) THEN
        CALL GETABT(PERLEN,NUMTS,TIMX,NPER,TIMREL,TBEGIN,TBGABS,KKPER,
     1            KKSTP,IERR)
 
      ELSE IF(KKPER.EQ.0 .OR. KKSTP.EQ.0) THEN
        TIMREL=TBGABS
        CALL GETPS(PERLEN,NUMTS,TIMX,NPER,TIMREL,TBEGIN,KKPER,KKSTP,
     1           IERR)
        IF(IERR.GT.0) THEN
          WRITE(*,*) 'SPECIFIED STARTING TIME IS OUT OF RANGE. STOP.'
          WRITE(I7,*) 'SPECIFIED STARTING TIME IS OUT OF RANGE. STOP.'
		call stopfile  ! emrl
          STOP
        END IF
C
      END IF
 
      WRITE(I7,*) ' '
      WRITE(I7,5310) KKPER,KKSTP,TIMREL,TBGABS
      WRITE(I7,*) ' '
5310  FORMAT(1X,'REFERENCE TIME DATA:'/5X,'STRESS PERIOD =',I4
     1/5X,'TIME STEP =',I4/5X,'RELATIVE TIME WITHIN TIME STEP = ',F7.5/
     25X,'ABSOLUTE REFERENCE TIME = ',E13.5)
      END IF
 
C
C... ASSEMBLE FLOW DATA. READ STRESS PACKAGE & BUDGET FILES. GENERATE
C    CBF IF NECESSARY
C
      CALL HQDATA(HEAD,QX,QY,QZ,BUFF,IUNIT,IBOUND,IBUFF,QSTO,QSS,DELR,
     1 DELC,
     2 LAYCON,KKPER,KKSTP,INHQ,FNAME,PERLEN,NUMTS,TIMX,IBSTRT,HDRY,
     3 HNOFLO)
C
 
C... INPUT PARTICLE STARTING LOCATIONS
C
C  READ PARTICLE LOCATIONS FROM DATA FILE
C
C... SCRATCH FILE TO HOLD INITIAL PARTICLE DATA & DISCHARGE CODES
      IRLI5= 24
      IF(IBYTES.NE.1) THEN
      IRM= MOD(IRLI5,IBYTES)
      IF(IRM.NE.0) THEN
      WRITE(I7,*) 'RECORD LENGTH NOT EQUAL TO EVEN MULTIPLE OF 2. STOP.'
		call stopfile  ! emrl
      STOP
      END IF
      IRLI5= IRLI5/IBYTES
      END IF
      OPEN(UNIT=I5,STATUS='SCRATCH',ACCESS='DIRECT',FORM='UNFORMATTED',
     1RECL=IRLI5)
C
 
 
C... OPTION WAS SELECTED TO READ STARTING LOCATIONS FROM A FILE.
C
      IZRO=0
      IF(IRP.EQ.1) THEN
      WRITE(I7,5000)
      WRITE(I7,*) 'STARTING LOCATIONS NOW BEING READ FROM FILE...'
      N=0
      ISP=0
      IPR=0
      TROLD= -1.0E+30
 
C... FIND OUT IF STARTING LOCATIONS ARE BEING READ FROM AN ENDPOINT FILE
      CALL SLFTYP(ISLTYP,I6,I7,IESLOP)
 
150   CONTINUE
      IF(ISLTYP.EQ.1) THEN
        READ(I6,'(A)',END=190) LINE2
        IF(LINE2.EQ.' ') GO TO 190
        ICOL=1
        CALL URWORD(LINE2,ICOL,IWSTRT,IWLAST,2,J,RDUMMY,I7,I6)
        CALL URWORD(LINE2,ICOL,IWSTRT,IWLAST,2,I,RDUMMY,I7,I6)
        CALL URWORD(LINE2,ICOL,IWSTRT,IWLAST,2,K,RDUMMY,I7,I6)
        CALL URWORD(LINE2,ICOL,IWSTRT,IWLAST,3,NDUMMY,XVALUE,I7,I6)
        CALL URWORD(LINE2,ICOL,IWSTRT,IWLAST,3,NDUMMY,YVALUE,I7,I6)
        CALL URWORD(LINE2,ICOL,IWSTRT,IWLAST,3,NDUMMY,ZVALUE,I7,I6)
        CALL URWORD(LINE2,ICOL,IWSTRT,IWLAST,2,JCODE,RDUMMY,I7,I6)
        IF(LINE2(IWSTRT:IWLAST).EQ.' ') THEN
          JCODE=0
          ICODE=0
          KCODE=0
          TRLEAS=0.0
        ELSE
          CALL URWORD(LINE2,ICOL,IWSTRT,IWLAST,2,ICODE,RDUMMY,I7,I6)
          CALL URWORD(LINE2,ICOL,IWSTRT,IWLAST,2,KCODE,RDUMMY,I7,I6)
          CALL URWORD(LINE2,ICOL,IWSTRT,IWLAST,3,NDUMMY,TRLEAS,I7,I6)
          IF(TRLEAS.NE.0.0 .AND. IREV.EQ.1) GO TO 150
          IF(TRLEAS.GE.TSTOP) GO TO 150
        END IF
      ELSE
        CALL EPSLOC(I6,I7,IESLOP,J,I,K,XVALUE,YVALUE,ZVALUE,JCODE,
     1              ICODE,KCODE,TRLEAS,IRQUIT,NCOL,NROW)
        IF(IRQUIT.EQ.1) GO TO 190
      END IF
C
C... IF RUN IS TRANSIENT, READ NEW HEADS IF NECESSARY
      IF(ISS.EQ.0) THEN
        IF(TRLEAS.NE.TROLD) THEN
          TRLABS= TBGABS + TRLEAS
          CALL GETPS(PERLEN,NUMTS,TIMX,NPER,TRLABS,TBEGIN,KPR,
     1               KSP,IERR)
          IF(IERR.GT.0) GO TO 150
          IF(KPR.NE.IPR .OR. KSP.NE.ISP) THEN
            ISP=KSP
            IPR=KPR
            CALL CBFHED(HEAD,NCOL,NROW,NLAY,I10,NSTEP,IERR,NUMTS,NPER,
     1                  NSTEPF,NSTEPL,KPR,KSP,2,I7)
          END IF
          TROLD=TRLEAS
        END IF
      END IF
C
      CALL GETIJK(I,J,K,XVALUE,YVALUE,ZVALUE,ZVL,ZV,ICODE,
     1 JCODE,KCODE,HNOFLO,XMAX,YMAX,ZTOP,ZBOT,HEAD,LAYCON,NCON,NCOL,
     2 NROW,NLAY,NZDIM,IGRID,HDRY,TRLEAS)
C
      IF(I.EQ.0 .OR. J.EQ.0 .OR. K.EQ.0) GO TO 150
C
      N=N+1
      INILOC(N)= (K-1)*NCOL*NROW + (I-1)*NCOL + J
      JLC(N)=J
      ILC(N)=I
      KLC(N)=K
      XLC(N)=XVALUE
      YLC(N)=YVALUE
      ZLC(N)=ZV
      ZLLC(N)=ZVL
      IF(TRLEAS.EQ.0.0) THEN
        TOT(N)= -1.0
        IDCODE=0
      ELSE
        TOT(N)= -1.0E+30
        IDCODE= -2
      END IF
C
C... STORE SOME PARTICLE DATA IN A DIRECT ACCESS SCRATCH FILE
      WRITE(I5,REC=N) XLC(N),YLC(N),ZLLC(N),IDCODE,TRLEAS,IZRO
C
      IF(N.LT.NPART) GO TO 150
190   CONTINUE
      NPRT=N
      WRITE(I7,*) ' STARTING LOCATIONS HAVE BEEN READ FROM FILE'
      IF(NPRT.EQ.0) THEN
        WRITE(I7,*) ' NUMBER OF VALID PARTICLES = 0. RUN STOPPED.'
		call stopfile  ! emrl
        STOP
      ELSE
        WRITE(I7,5300) NPRT
5300    FORMAT(' NUMBER OF VALID PARTICLES =',I10)
      END IF
 
      END IF
 
 
C... OPTION WAS SELECTED TO GENERATE STARTING LOCATIONS INTERACTIVELY
C
      IF(IRP.EQ.2) THEN
      MES= ' STARTING LOCATIONS WILL BE GENERATED FOR 3-D SUBREGIONS.'
      CALL MPRMPT(MES)
5180  FORMAT(1X,A)
      CALL MPRMPT(' ')
      MES= 'SHOULD PARTICLES BE PLACED IN CONSTANT HEAD CELLS ?'
      CALL YESNO (MES,ICHEAD,'2.1.46')
      N=0
      KO=0
      IREGEN=0
      IF(IREV.EQ.0) THEN
      CALL MPRMPT('SELECT AN OPTION FOR RELEASING PARTICLES:')
      CALL MPRMPT('  1 = SINGLE, INSTANTANEOUS RELEASE')
      CALL MPRMPT('  2 = MULTIPLE RELEASE TIMES')
      CALL GETI(1,IREGEN,3,1,2,'2.1.14A')
      IREGEN=IREGEN-1
      END IF
200   CONTINUE
      KO=KO+1
      WRITE (STRING,5190) KO
      CALL MPRMPT(STRING)
5190  FORMAT(' ENTER DATA FOR SUBREGION ',I2,' :')
      CALL MPRMPT(' ')
      CALL MPRMPT('DEFINE THE SUBREGION BY ENTERING THE MINIMUM AND MAXI
     1MUM J,I,K COORDINATES.')
      CALL MPRMPT('  ENTER:  MINIMUM J, MAXIMUM J, MINIMUM I, MAXIMUM I,
     1 MINIMUM K, MAXIMUM K')
      CALL GETI(6,IARRAY,0,0,0,'2.1.14')
      JMIN=IARRAY(1)
      JMAX=IARRAY(2)
      IMIN=IARRAY(3)
      IMAX=IARRAY(4)
      KMIN=IARRAY(5)
      KMAX=IARRAY(6)
      IF (JMAX.LT.1.OR.JMAX.GT.NCOL) JMAX=NCOL
      IF (IMAX.LT.1.OR.IMAX.GT.NROW) IMAX=NROW
      IF (KMAX.LT.1.OR.KMAX.GT.NLAY) KMAX=NLAY
      KWT=0
      IF(KMIN.EQ.0) THEN
      KMIN=1
      KMAX=1
      KWT=1
      END IF
      CALL MPRMPT('WHERE SHOULD THE PARTICLES BE LOCATED ?')
      CALL MPRMPT('  1 = WITHIN A CELL')
      CALL MPRMPT('  2 = ON ONE OR MORE OF THE CELL FACES')
      CALL GETI(1,ILOC,3,1,2,'2.1.15')
      ILOC=ILOC-1
      IF(ILOC.EQ.0) THEN
      CALL MPRMPT('ENTER A 3-D ARRAY OF PARTICLES: NX  NY  NZ')
      CALL GETI(3,IARRAY,0,0,0,'2.1.16')
      NJ=IARRAY(1)
      NI=IARRAY(2)
      NK=IARRAY(3)
      IF(NJ.EQ.0) NJ=1
      IF(NI.EQ.0) NI=1
      IF(NK.EQ.0) NK=1
 
      ELSE
      IF1=0
      IF2=0
      IF3=0
      IF4=0
      IF5=0
      IF6=0
      CALL MPRMPT('WHAT CELL FACES DO YOU WANT TO PLACE PARTICLES ON ?')
      CALL MPRMPT('  (ENTER ALL THE FACE NUMBERS ON A SINGLE LINE)')
      CALL GETSTR(CFACES,0,'2.1.47')
 
      ICOL=1
      LCFACE= LEN(CFACES)
300   CONTINUE
      CALL URWORD(CFACES,ICOL,IWSTRT,IWLAST,2,NFACE,RDUMMY,I7,IUNIT(1))
      IF(CFACES(IWSTRT:IWLAST).EQ.' ') GO TO 310
      IF(CFACES(LCFACE:LCFACE).EQ.'E') GO TO 300
      IF(NFACE.EQ.1) THEN
        IF1=1
      ELSE IF(NFACE.EQ.2) THEN
        IF2=1
      ELSE IF(NFACE.EQ.3) THEN
        IF3=1
      ELSE IF(NFACE.EQ.4) THEN
        IF4=1
      ELSE IF(NFACE.EQ.5) THEN
        IF5=1
      ELSE IF(NFACE.EQ.6) THEN
        IF6=1
      END IF
      GO TO 300
310   CONTINUE
 
      IF(IF1.EQ.1) THEN
      CALL MPRMPT('ENTER A 2-D ARRAY OF PARTICLES FOR FACE 1:  NY  NZ')
      CALL GETI(2,IARRAY,0,0,0,'2.1.17')
      NI1=IARRAY(1)
      NK1=IARRAY(2)
      END IF
 
      IF(IF2.EQ.1) THEN
      CALL MPRMPT('ENTER A 2-D ARRAY OF PARTICLES FOR FACE 2:  NY  NZ')
      CALL GETI(2,IARRAY,0,0,0,'2.1.18')
      NI2=IARRAY(1)
      NK2=IARRAY(2)
      END IF
 
      IF(IF3.EQ.1) THEN
      CALL MPRMPT('ENTER A 2-D ARRAY OF PARTICLES FOR FACE 3:  NX  NZ')
      CALL GETI(2,IARRAY,0,0,0,'2.1.19')
      NJ3=IARRAY(1)
      NK3=IARRAY(2)
      END IF
 
      IF(IF4.EQ.1) THEN
      CALL MPRMPT('ENTER A 2-D ARRAY OF PARTICLES FOR FACE 4:  NX  NZ')
      CALL GETI(2,IARRAY,0,0,0,'2.1.20')
      NJ4=IARRAY(1)
      NK4=IARRAY(2)
      END IF
 
      IF(IF5.EQ.1) THEN
      CALL MPRMPT('ENTER A 2-D ARRAY OF PARTICLES FOR FACE 5:  NX  NY')
      CALL GETI(2,IARRAY,0,0,0,'2.1.21')
      NJ5=IARRAY(1)
      NI5=IARRAY(2)
      END IF
 
      IF(IF6.EQ.1) THEN
      CALL MPRMPT('ENTER A 2-D ARRAY OF PARTICLES FOR FACE 6:  NX  NY')
      CALL GETI(2,IARRAY,0,0,0,'2.1.22')
      NJ6=IARRAY(1)
      NI6=IARRAY(2)
      END IF
      END IF
 
C
C... TURN OFF PARTICLE REGENERATION IF NECESSARY (OPTIONAL FOR FORWARD
C    TRACKING AND NOT ALLOWED FOR BACKWARD TRACKING)
      IF(IREGEN.EQ.0) THEN
        NRGEN=0
        DTR=0.0
        TIMINI=0.0
        TIMFIN=0.0
      ELSE
        CALL MPRMPT('SPECIFY THE TIME PERIOD FOR RELEASING PARTICLES:')
        CALL MPRMPT
     1('  INITIAL TIME, FINAL TIME, RELEASE INTERVAL, & TIME CONVERSION
     2FACTOR')
        CALL GETR(4,RARRAY,0,0.0,0.0,'2.1.53A')
        TIMINI=RARRAY(1)*RARRAY(4)
        TIMFIN=RARRAY(2)*RARRAY(4)
        DTR=RARRAY(3)*RARRAY(4)
C
        IF(TIMFIN.LT.TIMINI) THEN
          TIMINI=0.0
          TIMFIN=0.0
          DTR=0.0
          NRGEN=0
        ELSE
          IF(DTR.LE.0.0) THEN
            NRGEN=0
            DTR=0.0
          ELSE
            NRGEN= (TIMFIN-TIMINI)/DTR
            IF(NRGEN.LT.0) NRGEN=0
          END IF
        END IF
      END IF
C
      CALL GENPAR(N,NRGEN,TBGABS,TBEGIN,TIMINI,TIMFIN,DTR,ICHEAD,
     1  ILOC,IBATCH,ILSTOR,NI,NJ,NK,IF1,IF2,IF3,IF4,IF5,IF6,NI1,NK1,
     2  NI2,NK2,NJ3,NK3,NJ4,NK4,NJ5,NI5,NJ6,NI6,NPART,NCOL,NROW,NLAY,
     3  ISS,KWT,IBOUND,HEAD,JLC,ILC,KLC,XLC,YLC,ZLLC,HNOFLO,HDRY,
     4  PERLEN,NUMTS,TIMX,TOT,INILOC,XMAX,YMAX,NPER,I5,I6,I10,I7,
     5  JMIN,JMAX,IMIN,IMAX,KMIN,KMAX,TSTOP)
C
      WRITE (STRING,5200) N,NPART
      CALL MPRMPT(STRING)
5200  FORMAT(1X,I6,' PARTICLES USED OUT OF ',I6,' MAXIMUM')
      MES= 'DO YOU WANT TO SPECIFY ANOTHER SUBREGION ?'
      CALL YESNO(MES,IANS,'2.1.53')
      IF(IANS.EQ.1) GO TO 200
      NPRT=N
      IF(ILSTOR.EQ.1) CLOSE (I6)
      END IF
C
C  CHANGE ZONE CODES IN THE IBOUND ARRAY
C
      IF (IZSTOP.NE.0) THEN
      MES= 'DO YOU WANT TO CHANGE ANY OF THE ZONE CODES IN THE IBOUND AR
     1RAY ?'
      CALL YESNO (MES,IANS,'2.1.54')
      IF (IANS.EQ.1) THEN
250   CONTINUE
      CALL MPRMPT('WHAT TYPE OF CHANGE DO YOU WANT TO MAKE ?')
      CALL MPRMPT('   1 = CHANGE AN ENTIRE LAYER')
      CALL MPRMPT('   2 = CHANGE AN INDIVIDUAL CELL')
      CALL MPRMPT('   3 = CHANGE ALL CELLS IN A SUBREGION OF THE GRID')
      CALL GETI(1,IANS,3,1,3,'2.1.23')
      IF(IANS.EQ.1) THEN
      CALL MPRMPT('ENTER THE LAYER NUMBER:')
      CALL GETI(1,NLAYER,3,1,NLAY,'2.1.24')
      CALL MPRMPT('ENTER THE NEW ZONE CODE:')
      CALL GETI(1,NUM,1,1,0,'2.1.25')
      DO 260 I=1,NROW
      DO 260 J=1,NCOL
      IF(IBSTRT(J,I,NLAYER).NE.0) IBSTRT(J,I,NLAYER)=NUM
260   CONTINUE
      ELSE IF (IANS.EQ.2) THEN
      CALL MPRMPT('ENTER THE CELL INDICES:  J  I  K')
      CALL GETI(3,IARRAY,0,0,0,'2.1.26')
      J=IARRAY(1)
      I=IARRAY(2)
      K=IARRAY(3)
      CALL MPRMPT('ENTER THE NEW ZONE CODE:')
      CALL GETI(1,NUM,1,1,0,'2.1.27')
      IF(IBSTRT(J,I,K).NE.0) IBSTRT(J,I,K)=NUM
      ELSE IF (IANS.EQ.3) THEN
      CALL MPRMPT('DEFINE THE SUBREGION BY ENTERING THE MINIMUM AND MAXI
     1MUM J,I,K COORDINATES.')
      CALL MPRMPT('  ENTER:  MINIMUM J, MAXIMUM J, MINIMUM I, MAXIMUM I,
     1 MINIMUM K, MAXIMUM K')
      CALL GETI(6,IARRAY,0,0,0,'2.1.28')
      JJ1=IARRAY(1)
      JJ2=IARRAY(2)
      II1=IARRAY(3)
      II2=IARRAY(4)
      KK1=IARRAY(5)
      KK2=IARRAY(6)
      IF(JJ2.LT.1.OR.JJ2.GT.NCOL) JJ2=NCOL
      IF(II2.LT.1.OR.II2.GT.NROW) II2=NROW
      IF(KK2.LT.1.OR.KK2.GT.NLAY) KK2=NLAY
      CALL MPRMPT('ENTER THE NEW ZONE CODE:')
      CALL GETI(1,NUM,1,1,0,'2.1.29')
      DO 270 K=KK1,KK2
      DO 270 I=II1,II2
      DO 270 J=JJ1,JJ2
      IF(IBSTRT(J,I,K).NE.0) IBSTRT(J,I,K)=NUM
270   CONTINUE
      END IF
      MES= 'DO YOU WANT TO CHANGE SOME MORE ZONE CODES ?'
      CALL YESNO (MES,IANS,'2.1.55')
      IF (IANS.EQ.1) GO TO 250
      END IF
      END IF
C
      RETURN
280   FORMAT(A)
      END
